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Abstract 

A new analytical operator method is discussed which solves linear ordinary differential equations with regular 
singularities. Solutions are obtained in analytic series form and also in Mellin-Barnes-type contour integral 
form. Exact series solution is obtained without having to calculate series coefficients by recurrence relation. Both 
homogeneous and inhomogeneous equations are solved identically without having to calculate the Green's function 
explicitly in the case of inhomogeneous equation. Closed-form solutions are obtained for all the special functions 
appearing in mathematical physics. For a second-order equation both the independent solutions are obtained 
without invoking Wronskians, even when the indices differ by an integer. 

Keywords: special functions, exact solutions, contour integral solution, analytic series solution 
1 Introduction 

Operator methods for obtaining solutions of ordinary differential equations with constant coefficients are well 
known and well studied. However, there are very few operator methods in the literature to solve linear ordinary 
differential equations of all orders with variable coefficients. The normal method to solve homogeneous equations 
with regular singularities is the Frobenius Method, where we assume a series and put it in the equation to 
obtain a recurrence relation between the series coefficients. However we face difficulties when the indices at 
the regular singularities differ by an integer. The method involving Wronskians generally applied to such cases 
is quite complicated. For inhomogeneous equations, calculation of Green's function using the solutions of the 
corresponding homogeneous equation is also difficult. There are really very few methods which treat homogeneous 
and inhomogeneous equations in a simple and unified way. 

To provide such a unified method to solve ordinary and partial, homogeneous and inhomogeneous, linear and 
nonlinear differential equations, the Adomian decomposition method PQ was designed. Using the inverse of 
the highest derivative, the differential equation is converted into an integro-differential equation, which is then 
solved by a successive approximation-type method to give an analytic series solution. While this method works 
beautifully in many cases, it fails severely when regular singularities are encountered. Out of a few thousands 
of papers published on the convergence and verification and modification of Adomian's method, very few really 
acknowledge this fact. One such is [2], where a modification of the invertible operator is proposed which takes 
care of the singularities. However, the authors have not given any particular way to solve the problem and have 
suggested that the form of the equation should decide the method. As an example they show how a lower-order 
operator (contrary to the inverse of the highest derivative in Adomian's method) can be used to obtain Bessel 
polynomials from the defining second-order differential equation, which has an irregular singularity. Though 
the authors have been successful in obtaining one solution , there is no way to obtain the other one using their 
method. Also, no justification has been provided for using Picard's iteration on an integro-differential equation 
with a singular kernel. 



In this paper we shall give a general method to obtain solutions of ordinary linear differential equations having 
regular singularities. Using a particular kind of operator the differential equation shall be converted into an 
inhomogeneous integro-differential equation, the inhomogeneous part being just a function of initial conditions 
and the original inhomogeneous function, if present from before. The resolvent of the integro-differential operator 
acting on the inhomogeneous part gives the solution. 

The integro-differential operator shall be shown to be bounded and contractive so that the well known Banach 
theory can be used to obtain a binomial series expansion of the resolvent operator. The infinite series of the 
operator then acts on the inhomogeneous part to generate the series solution. 

In the next section we shall use a contour integral representation of the resolvent operator to obtain Mellin- 
Barnes-type contour integral solutions to differential equations. 

2 Operator method to solve linear ordinary differential equations 

The operator method can be used to solve general nth-order linear differential equations. We shall work mainly 
with second-order ODEs in this paper and cover the most important special functions of mathematical physics 
which satisfy second-order ODEs. However the methods can be easily generalised to solve higher-order ODEs. 
We first study equations with a regular singularity at the origin and another (regular or irregular) singularity at 
infinity. 

2.1 Two singularities at and oo 

The general second-order equation with singularities at and oo is: 

d 2 ^{z) djj(z) 
dz 2 P ( Z '~dz~ + q ( z >^ z > = F W- W 

If p(z) and q(z) have series expansions p(z) = Y17L~iPi zl anc ^ l( z ) = J2iL-2 1i z% i then the singularity at z = is 
regular [12] . We take a transformation ip(z) — z x f(z). Such a transformation enables us to remove the term z~ 2 
from the coefficient of f(z). Equating the coefficient of z~ 2 to zero, we get the indicial equation 

A 2 + (p_i-l)A + $_ 2 = 0. (2) 

Let Ai, A2 be the roots of the indicial equation, with Ai > A2. If Ai — A2 is an integer then the second solution 
would have a branch-point singularity or logarithmic singularity. The equation for a particular index A, satisfied 
by f(z), is 



^ + (f + E + (E Di**)M = *- X F{z), (3) 

8=0 ' i=0 

where a = 2A +p_i,Cj = pi and Di = Xpi + qi—i- The above equation can be written as 



hi^^h + «E^)^r + £>*')/(*)) = (4) 

i=0 i=0 

We define a new operator 

L = J ^(J z a (-)dz)dz, (5) 
where the integrals are indefinite. Operating L on (U), we get 
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/(*) + / ^{j z\{j^C i z i )^ + {Y j D i z i )f{z))dz)dz = J ±(J ' z a - x F(z)dz)dz + h{z). (6) 



i=0 i=0 



fo(z) is the 'integration constant' of the operator L, satisfying the constraint 



i_d 

z a dz y dz ' w 

From (|Tj) , we get 

/o(*0 = c + ci y ~^ dz - (8) 

We now define another operator 

i= i=0 

Thus, (JSJ) can be written in operator notation as 

(f + A)- f(z) = L-z- x F(z) + f (z). (fO) 
The solution f{z) of (flU)) can be formally written as 

/(z) = (1 + A)- 1 • L ■ z- x F{z) + (1 + A)- 1 • f (z). (11) 

From the above equation, we can see that first term on the right-hand side is the particular solution while the 
second part is the complementary solution, which is also the solution of the corresponding homogeneous equation. 
The Green's function operator can be seen to be (1 + A)^ 1 ■ L ■ z~ x Q. 

2.2 Three singularities at 0, 1 and oo 

The general second-order equation with singularities at 0, 1 and oo can be written as: 

z) fm +p{z) wp_ + q{mz) = F{zy (12) 

az z dz 

If p(z) and q(z) have series expansions p(z) = z X^-i Vi z% an d = z ^^t = -2 ( l iZ% i f nen we have a regular 
singularity at 0. If the differential equation has singularities at any 3 other points, then the singularities can be 
mapped onto 0, 1 and oo without changing the indices, as explained in [12] . 

As before, we take a transformation ip(z) = z x f{z) to remove the term z~ 2 from the coefficient of f(z). We 
obtain the indices by solving the indicial equation likewise. The equation for a particular index A, satisfied by 
f(z), is 

—^ + (- + l^ C ^)— + C^D l z)f{z)-z- [ ^ = z F(z), (13) 

i=0 i=0 
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where a = 2A +p~%, C — p - 2A, D = A(l - A) + Xp + Z?i = A(l — A) + Api + go and A = Ap^ + <&_i, 
Ci = pi for all other i. The above equation can be written as 



hi <^> + «£ «^ + E - ^) - ^.). (») 

i=0 i=0 

Defining L as in (O and operating L on (fT4)l . we get, as before, 

(l + A)./(z) = L.z-V(z) + / (z), (15) 

where /o(-z) is as defined before in ([8|), and 

A= f -!( f z^-z^ + ^CiZ^ + ijrD^^dz. (16) 

J Z J Z i—o Z i=0 

The solution /(z) can again be written as (fTTj) 

In the next section, we shall study the properties of the operators L and A. 
3 The operators L and A 

Let G be an open subset of C defined by {z S G : < \z\ < 1}. Let v : G — > R+ be a continuous and strictly 
positive function on G, here called the weight. We define the following weighted Banach space of holomorphic 
functions on G: 

Hv(G) := {/ e H(G) : ||/||„ := sup v(z)\f(z)\ < oo} (17) 



We refer to [5J O El IHj , and the references given therein, for details on the above Banach space. From the standard 
literature on solutions of second-order ODEs, it is obvious that our space contains functions of the regular and 
irregular kinds. The regular kind is finite at the origin and has a Taylor series expansion about the origin: 

oo 

Mz)=^2cr i z i . (18) 

i=0 



The irregular one has multiple pole and even branch-point singularity at the origin 

'z- 5x EZo cn ^ 1 if(Ai-A 2 )gN, 
f ir (z) = { (logz) if (Ax = A 2 ), (19) 



- 5X EZo cn ^ + ( lo § z ) TZo if ( A i - A 2) g n, 



where 5 A = Ai — A 2 = a — 1. 
We choose v(z) — \z\ a , so that 



lim f(z)v(z) = (20) 



for both regular and irregular functions. 

Theorem 3.1. There is a subspace Gq = {z G Go : < \z\ < zq\ C G such that 

i. L : Hv(Gq) — > Hv(Gq) is continuous, 

ii. A : Hv(Gq) — > Hv(Gq) is continuous, and 

in. (1 + ^l)- 1 = 1 - A + A 2 - A 3 + A 4 + ... + (-l)J'AJ' + ... 
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Proof, i. For a^Nwe have f(z) = J2m=sx c mZ m - 

From the definition of L in ©, we have: Lz m = J ^(/ z a+m dz)dz = (m+2 )^ + + 2 „ t+1) ■ Therefore ||£« m || v < 
M|b m |L where M = ? -nm^r < 1- 

For a € N we have integral powers of z starting from z 1 ~ a (for a ^ 1) as well as z" l log(2;) terms in f(z) with 
to = 0, 1, 2..oo. For a = 1, 

L- z m = (m +2 )2 and Lz m log(z) = ^fp- - 2 (m+2)a • 
For a £ N but 1 and to — 0, 1, 2..oo, 

i i f? — i nwo i v if m ^ a — 3, 

Lz = < ( m + 2 )( 3 + m ~ a ) 

) og ^s if m = a — 3. 



For the z™log(z) terms with n £ K+, 

i-z m log(z) = z^ +m ( 



. («-l) 



2+ -3 -2m -a ( log(z) 



(2 + m) 2 (l + m + a) 2 (2 + m)(l + m + a) ' 



Thus, from all the above cases, we can see that for / defined as in ([TBI fT9| . L : / — > /. Also, ||L/(z)||„ < ||/(z)||d 
for all < |z| < |zo| < 1. This proves that L is continuous. 

ii. For A defined as in ©, we can immediately see that for a £ N, 

A . Z m = V ( m + A )P' „i+m+l 

(i + to + a)(i + to + 1) 



z 



The series obtained is seen to be convergent from the ratio test. Therefore, ||Az m ||,j < M||z m ||.„, where M 

(m+X)pi- 

For the z n log(z) terms with n £ 



sup£~o /(*) = E™ = -,a CmZ™, therefore ||A/(z)||„ < M||/(2)||„. 



T v m \r*a-( -A — -(Api+9i-i)(2i+2fc+a+l)+a+i--fc +a(l+i) j+ m +l , y^oo (m+A)p i +g i _ 1 ,j+ m +li„„^\ 

Zv • Z tO&\Z) — l^ l=0 (i+m+a)*(i+m+iy 2 Z " t " Z^i=0 (i+m+a)(i+m+l) Z iu 6l 2 J- 

For A defined as in (fT6|) . we have just another extra term L(z d ^ 



dz 2 r 

Uz^f^) = . ^ m+1 fOT ™ ^ -1 and to > 1 - a 

v. dz 2 / (m+l)(m+a) ' — 

For to = -1, L(z^pl) = ^]o g {z). 

T („ d 2 (z m log(z))) , _ (-a+2m tt +m 2 (2+a)+m(- 1+m 2 ) (m+a) log( Z )) m+1 f , , 

M Z 3? ) — (m+l)2(m+a)2 Z ' 

From the above, it is seen that A: f —> f and H-L/^)^ < M||/(z)|| t ,. 

iii. We can choose zq such that M < 1. Therefore, for 0<|z|<|zo|<l,^4isa contraction operator and hence 
we can expand the resolvent (1 + A)" 1 in a Neumann series. 

□ 
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4 Representations of the resolvent operator 

In the last section, we saw that the resolvent operator (1 + A)^ 1 can be expanded in a Neumann series for 
< \z\ < \zo\ < 1, so that M < 1. We can also have a contour-integral representation of (1 + A)^ 1 : 

ra-\-LOO 

2m _ 

where < a < 1 and x{z) is an element of the Banach space Hv(G) 



(1 + A)~ 1 x(z) = ^- I r(s)r(l-s)(-A)-*x(z)ds, (21) 



This form of the resolvent operator can be obtained by taking an inverse Mellin transformation of the expression 
for the fractional power of a closed operator given by V. Balakrishnan in [5] and [TT] . 

We consider the expression $(t) = ^7 /(? r(s)r(l ~ s)(—t)~ s ds, where C is a Mellin-Barnes-type contour. For 
\t\ < 1 we close the contour on the left side of the complex plane, and so the singularities enclosed by the contour 
are those of T(s) at s = —n, n E Z+. Using Cauchy's residue theorem, we obtain $(<) = Y^Lo^ 1 - Similarly, the 
expression for the resolvent operator defined by (|21[) gives the Neumann series for \ \A\\ V < 1. 

In the next section, we shall apply the operator method to obtain exact solutions of various well-known differential 
equations occurring in mathematical physics. 

5 Examples from Mathematical Physics 
A. Homogeneous equations 
5.1 Singularities at 00 

To begin with, we solve the most elementary differential equations with constant coefficients. 
5.1.1 The exponential function 

The first-order ODE satisfied is -j| = y(z).We integrate both sides of the equation to get 

(1 - J)y(z) = c, (22) 

where J is the indefinite integration operator and c is the integration constant. Thus, our operator A is 3. For 
\z\ < 1, < 1. The solution is 

y(z) = (l-3)- 1 c. (23) 
We first expand the resolvent in a Neumann series. So, 

2 3 

y(z) = (l + 3 + 3 2 + J 3 ...)c = c{l + z+^+ Z — + ...) = ce z . (24) 

Next we use the contour- integral form of the resolvent. A fractional power of 3 is identical to that of the definite 
integral with limits and z. We know from fractional calculus that 

Therefore, 

1 />a+too 1 />a+too 

(1 - 1)- 1 ? = 5- / r(*)r(i - s)(-3y s cd s = —c / r( s )(-z)- s ds. (26) 

The last expression is the inverse Mellin transform of gamma function, which is the exponential function. 
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5.1.2 Trigonometric and hyperbolic functions 



The second-order differential equation satisfied by these are ^| ± lu 2 z = 0. We choose L = J 2 . Acting both sides 
by L, we get 

(l±3 2 u 2 )y(z) = c + c 1 z, (27) 

where c®andc\ are constants. 

Therefore A is ±3 2 uj 2 . We choose \z\ < u^ 1 to ensure that < 1. The Neumann scries expansion method 

proceeds exactly as before. We get cos(z) and cosh(z) with c , and sin(z) and sinh(z) with c\. 

We have 

JZ - C T(l + q+P )- (28) 
Substituting p = —2s, q = for the term with cq, and q = 1 for the term with ci, we get 



y(z) = (1 + Ay^co + az) 

= — / r(*)r(i - s)(tu 2 3 2 )- s (c + c lZ )ds 

J a - LOO 

The above contour integrals, though not the exact inverse Mellin transforms, can be shown to yield the same 
trigonometric / hyperbolic functions, using Cauchy's residue theorem. 

5.2 Singularities at and 00 

We treat the most important special functions having singularities at and 00 — the Bessel functions and the 
confluent hypergeometric function. 

5.2.1 The Bessel functions 

This equation has a regular singularity at z — and an irregular singularity at 00. We take tp(z) = z K f{z). The 
indices at the regular singularity are k — ±v. The equation for f(z) for the positive index is 

which can be written as 

' V^ 1 ^) = -/(*)• (32) 



z 2l/+1 dz K dz 
We choose the operator L = J z2 l+i (J z 2v+1 ()dz)dz. 

The 'integration constant' fo(z) for L satisfies 



1 d 2v+1 djp_ 

i 2 »+ 1 dz dz ' x ' 
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A simple integration of the above equation gives fo(z) = c + c\ J z £ z +i . 

Acting on both sides by L, we get (1 + A)f(z) = fo(z), where A = L is a bounded operator. We choose z$ suitably, 
so that V \z\ < \z \, \\A\\ V < 1. 

The regular solution is obtained by acting the resolvent operator on c . Wc start with the Neumann series as 
usual. 

A ■ c = c 0T ^ (f ) 2 , A 2 ■ c = c 2(1+ „) (2+i/) (f ) 4 -. 
It can be proved by induction that 

An „ 1 (" z\2n 

Therefore, the regular solution is 

fr(z) = (1 - A + A 2 + ...)co 

which is just c T(l + v){^)~ v J n (z). 

Since the calculation of fractional powers of operators is difficult, a heuristic method to obtain the fractional power 
shall be followed from now on. The method is to calculate A n for n £ N and then to generalise the expression so 
that n can be fractional too. 



In this case, we see that writing 

A n • c = c 



1 I Z 



n\{l + v) n \2 



- C °T{l + v + n)T{l + n)\2) (35) 
solves our problem, as the T{n) functions are valid for both integral and nonintegral n. Therefore, 



fr(z) = — / r( s )r(i - S )(-A)- S c d s 

(l I ^ ds. (36) 



—c f 



When v ^ R+, then the second irregular solution can be easily obtained by acting the resolvent on the term with 

Ci. 

A . r , z~ 2 " — _£3_^L_ (z\2 A2 . z~ 2 " _ _ci 1 I z\i 

A Cl -2v ~ -2v 1-v \2> ' 71 Ll -2v ~ -2v2(\-v)(2-v) \2> — 

For general v, 



2v-2v 



1 ' ) = (-1) " r(i-!- + «)r(i + f) (,iJ 
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Thus we obtain the irregular function as 



2„ ci / 1 / z\ 1 /r 



-2i/ V V 2 / 2(1 — i/)(2 — f) V2 

1 Cl r-,„ r (s) /_*y 2 ^ (38) 



27rt -2i/ y a _ loo r(i-i/- s )V 2 

Which is ^r(l - !/)(§)-" J_ n (z). 

In fact equations (|35l37p are special cases of 



(-l) v x« fx"" 2 



^> = (i +f ). ( i +f+ »). UJ (39) 

where (a)„ is the Pochhammer symbol. 

When v = n € R+, logarithmic terms occur. We need a result from fractional calculus: 

T (z A ln(z)) = + ^ z^+"(ln(z) + ^(A + 1) - yj(X + v + 1)), (40) 

1 (A + v + 1) 

where ip is the digamma function with tp(l) — —7, the Euler constant. 

We need to obtain A~ s J z 1+2,1 . We observe that equation (|37|) is valid for v = n with v < n — 1. 
For u = n, we have 

.«/ f 1+3™ _ {-IT Hz) Ul) 

A{ J Z >-(2^ + nW (41) 

A n+1 (J z 1+2n ) = (-l)"-i_(|) 2 (]ii(.?) - J^r), A m +"(JV+ 2n ) = (-l)"(|) 2m (ci(m) ln(z) + c 2 (m)) for integral 
m. 

Therefore, 



A m+n+ll I z l+2n\ = (_!)"(£ 



2(m+l) 

2 

ci(m) , , . 2(1 + m)(l + m + n)c 2 (m) — cUm)(2 + 2m + n) \ 

■ ln ( z ) + on I Z3a7i r~ZTT~Ivl ( 42 ) 



(1 + m)(l +m + n) 2(1 + m) 2 (l + m + nf 



From the above, we obtain 



ci(m) 

cim + 1 =- £—t -, 43 

v ; (1 + m)(l + m + ri) 

, . 2(1 + m)(l + m + n)c 2 (m) - ci(m)(2 + 2m + n) , AA . 

c 2 (m+l = — — . 44 

2(1 + mj^(l +m + 



We uniquely solve (|43I44|) with the condition (|4Tj) to get 
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Cl{m) ( 4) "(l + n)r(l + m)r(l + n) 2 r(l + m + n)' (45) 

, )= (-l)T(2 + n) 

C2(m) 2( 1 + 2 ")(l + n)r(l + m)r(l+7i) 2 r(l + m + 7i) 1 ' 
(7 + + m) — ip(l + n) + ip{l + m + n)), 

A m+n ( f i+2n\ _ (~l) n (2 m(z) - 1 + j,(n + I) - ijj(m + 1) - ^(m + n+ 1)) /z\ 2m 
V7 / 2( 1 + 2 ™)r(n+ l)r(m+ l)r(n + m+ 1) V 2 / ' 



Using (|47p and putting v = m + n, we get 



/ / z i+2„\ = (-1)" (2 H*) - 7 + V»(n + 1) - ^(v - n + 1) - j,(v + 1)) /z\ 2{v ~ n) 

\J ) 2( 1 + 2 ")r(n+l)r(«-n+l)r(i; + l) ' 1 ' 

This equation is valid for fractional v as well. The irregular solution for integral n can thus be written as 

/r(z, =^ir rwr(i - s)Mr *(/^) d * 

1 21n(z) - 7 + V^(n + 1) r +t °° r ( s ) /zr 2(s+ "', 

as 



2m r(n+l) 7 a _ too 2 2 «+ir(l- s -n) V2 

1 1 /• a+, ' 00 r(s)(i/;(l- l s-n)+^(l-s))/^ x_2(s+ " ) 



2^r(n+l)7 Q _ too 2 2 «+ 1 r(l-s-n) 



as (49) 



5.2.2 The confluent hypergeometric function 



dz \ z I dz z 



The indices are and 1 — c. The ODE can be written as 



i*W^_ffi + !!/w|=0i (51) 



: rfz V dz J \ dz z 

We choose the operator L = J js(J z c ()dz)dz. The 'integration constant' fo(z) for L, satisfying 4- (z c ) = 
0, is fo(z) =c Q +J z- c dz. Operating by L, we get (1 - A)f(z) = f (z), where A = J £(/ z c ((^ + f {)))dz)dz. 

Lemma 5.1. For q e M+ but q N, 

(q + a) m _ r(a + m + g)r(l + q)r(g + c) g+m 



A m (z«) = ' , z 



(q + c) m (q + l) m r(a + q)T(c + m + q)T(l + m + q) 



z* 



Proof A(z q ) - (g+ a ) z<i+1 A 2 ■ cn - cn (■z+ a )('?+ a + 1 ) ~q+2 



(«+c)(l+<z)> " " u ^ ° u 2(g+ C )(g+c+l)- 

We assume A m (z q ) = ( g +e) + (]+i) z g+m , where (a)„ is the Pochhammer symbol. Acting v4 on it, we get 
A m+1 (z q ) = ^q^y~~fe^x) +1 z g+m+1 . Hence, the lemma is proved by induction. □ 



10 



From lemma [STTl we see that A is bounded and is a contraction operator for \z\ < 1. Putting g = in lemma [5TT1 
and expanding the resolvent (1 — A)~ x in a Neumann series, we obtain the regular solution : 

(a) re / £^ 
(c)„ V n! 



/r(z) - (1 + A + yl 2 + ...)c = c £ — = iFxfa, c, z). (52) 



The contour integral representation is immediately obtained as the expression for A n is valid even when n is 
fractional or complex. So, 

-i />a' + too -i -p/ \ />a' + (.oc -p/ \"p/ \ 

/»•(*) - IT- / T( s )r(l - s)(-A)-°cods = ^-^co / y la ' SJ (-z)-^, (53) 

27rt -/a'- t oc 27Ttr(a) ,/ a ,_ too r(c-s) 

with a' € (0,1). 

For c ^ R+, J z~ c dz — and therefore, for the irregular solution, we can still use lemma ISTTI with q = 1 — c. 

We Thus obtain 

fi(z) = (£ A n )z l ~ c = z 1 -* £ ( f^i) = ^"Wl + a - c, 2 - c, z). (54) 

The contour integral representation is 

Mz) . f ' + '°° r W r(i + a -c-.) 

2tti r(l +a-c) J a ,- LOO T{2-c-s) 

5.3 Singularities at 0, 1 and oo 
5.3.1 Gauss's hypergeometric function 

(1 - Z )^fi +(--(a + b + l)) d -^l- %{z) = 0. (56) 
dz \z J az z 

The indices at are and 1 — c. The ODE can be written as 

i dm f^m\ _( {a+b+1) m + ± f(z) + = o _ 



z 



z c dz \ dz J \ dz z dz 

We choose the operator L = f ■%(/ z c ()dz)dz. The 'integration constant' /o(^) for L, satisfying ~? j- {z c ^di^ ) 
0. is /o( z ) = c o + / z~ c dz. Operating by L, we get (1 — yl)/(z) = fo{z), where 

A = f±J*>(((a + b + l)§ + a l + z^)dz)dz 

Here onwards, the analysis is exactly similar to the previous (confluent hypergeometric) case. We just give the 
final results. 

A m ( q) = (q + a)m(q + b) m q+m = T(a + ro + q)T(b + m + q)T(q)T(q + c) q+m 

[ ' (q + c) m (q + l) m T(a + q)T(b + q)T(c + m + q)T(l+m + q) ' 1 j 

The regular solution obtained using (|55|) with g = 1 — c in the Neumann expansion of (1 — A)^ 1 is 

/r(*) = (f; A")c = co £ (^) = 2 F 1 (a, b, c, z). (59) 

1=0 n—1 ^ ' V * / 



n 



The contour integral representation is 



fr{z) = — / T( s )r(l - S )(~A)- S c ds 

27Ti J a' -loo 

i r(c) r( 5 )r( a - s )r(6- s ) 

= 7^— c , , , , / —, r -z) ds, (60 

27ri r(a)r(6) 7 Q ,_ too r(c-s) 



with a' € (0, 1) 



For — c ^ N, J z c = an( i therefore the irregular solution obtained using with q = 1 — c in the Neumann 
expansion of (1 — A)^ 1 is 

(l + a-c) n (l + b-c) n fz n 

i=0 n=l 



/r(z) = £ ^")co = co J2 (2 - c )„ ^ ) = 2 Fx (1 + a - c, 1 + b - c, 2 - c, Z ). (61) 



The contour integral representation is 

w „ = r '^ +t w + .-.-.ro +t -.-.) 

2tti r(l + a-c)r(l + 6-c) 7 a ,_ too T(2 - c - s) 

5.3.2 The Gegenbauer polynomials 

The differential equation satisfied by the Gegenbauer polynomials, which are general cases of Legendre polyno- 
mials, is 



(z 2 - 1)^-Pfi- + 2(0 + l)z^ - a(a + 2/3+ l)f(z) = 0. (63) 
dz z az 



A transformation z — > 2z — 1 is needed to bring the singularity at z = — 1 to z = 0. 
The transformed equation is 



:(l -z)^^ + G9 + l)(2z-l)^-,,m+2;-ll/(-.i=ii. ((,4! 



The solution to this equation follows exactly as in the case of Gauss's hypergeometric function. 

B. Nonhomogeneous equations 

5.4 Inhomogeneous Bessel equation 

2 d 2 4>{z) dMz) ,„ zA . 2 1 -^z y + 1 

*^ + + (1 - ^) = VHXIT^) (65) 

This inhomogeneous ODE has a particular solution which is the Struve function H u (z). Using equation (|11[) we 
calculate the particular solution to this ODE. 
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WW = z v (l + A)- 1 ■ L ■ ( ( 66 ) 



Using the fromulae from section (|5.2.1|) we obtain 



Using d2U) with q = 1 and (jST)) 



T^TTT J = ( ^ r ,3^ J (67) 



2 l ~ v z- 1 ., 2~ v fz^ 2v+1 



^'^f(lT^»°'-"' r(| + , + »)m + ,) UI < 68 > 



The series solution and contour integral solution of ip p {z) are therefore: 



= H v (z) 



-1)" 



2v+l+v 



1 

2ttI 



a' + too 



a — loc 



(-i) s r( 3 )r(i-.s) 
r(| + ^-,s)r(|- s ) 



1+I/-2S 



ds. 



(69) 
(70) 



Thus we can see that complicated Green's function method to solve inhomogcncous method can be avoided using 
this method. 

6 Conclusion 

In this paper a new operator method is studied which is capable of giving not only analytical series solution but 
also exact contour integral solution to all the special functions. Calculation of complicated Green's functions is 
avoided in this method while solving inhomogeneous equations. Though there is a restriction \z\ < 1 mentioned 
in this paper, we know that a Taylor series expansion remains valid till the nearest singularity of the differential 
equation and hence the solution obtained can be used well outside \z\ = 1 if the singularity is at infinity. Only 
the most important equations occuring in Mathematical Physics, with 2 or atmost 3 singularities. However this 
method can easily solve equations with more than 3 singularity like Heun's equation much in the same manner. 

In this paper all the differential equations considered have polynomial coefficients so a operator of the form [5] is 
used. In another paper this method shall be used to solve periodic differential equations with L = d ^ja , where 
D is the normal differentiation operator and lo, a frequency parameter. 

Only analytical solutions are discussed in this paper. In another paper we shall show the numerical efficiency of 
the method to solve eigenvalue problems not only of 2nd order but also higher order equations. This method has 
also been recently generalised to solve ODEs with algebraic non linearities. The operator in that case takes the 
form of an infinite operator continued fraction. 
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